#include <iostream>
#include <fstream>
#include <gsl/gsl_math.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <sstream>

#include <sys/stat.h>

#include "HeatBath.h"

//#include "IOParams.h"
#define DIM 2
#ifndef GAS_H_
#define GAS_H_

using namespace std;

int indx(int , int );
class Gas{
	public:
		int N;
	  double *theta;
		double *omega;
		double *vec_m;
		double *vec_M;
		double *force;
		double Ek;
		double Ep;
		
		Gas(int const);
    ~Gas();
    
    void initCoord();
    void initVeloc();
    
    void initCoord2();
    void initVeloc2(double Temp);
    
    void initCoord3();
    void initCondFromFile(string filename);
    
    // Calculate functions
    void calculateVec_m();
    void calculateVec_M();
    void calculateForce(double eps);
    
    // Heat reservatories.
    void addHeatBathTo(int partID, double Temp, double dt, double gamma, gsl_rng *r);
    
    // Energy
    void calculateEkin();
    void calculateEpot(double eps);
    
    // Update functions
    void updateCoord(Gas*, double dt, double dt2);
    void updateCoordHBLangevinTo(int part_id, Gas* oldgas, double dt, double dt2, HeatBath* hb);
    void updateCoordWhitoutHB(Gas* oldgas, double dt, double dt2);
    void removeCyclesCoord();
    void updateVeloc(Gas*, double dt);
    void updateVelocHBLangevinTo(int part_id, Gas* oldgas, double dt, HeatBath* hb);
    void updateVelocWhitoutHB(Gas* oldgas, double dt);
    

		// Save data to file
		void writeMacroState(ofstream *ptrfile, int ttime, double dt);
		void writeMicroState(ofstream *ptrfile, int ttime, double dt);
    
    // print functions
    void printCoord();
    void printVeloc();
    void printVec_m();
    void printVec_M();
    void printForce();
};

#endif /* GAS_H_ */
